Here we will recreate the ZooMSS model (version 2) in Heneghan et al. (2020 in preparation) using mizer.
We begin with some setup of required packages.
#get required packages
library(devtools)
## Warning: package 'devtools' was built under R version 3.6.3
## Loading required package: usethis
## Warning: package 'usethis' was built under R version 3.6.3
#most up to date master branch of mizer
#install_github("sizespectrum/mizer")
#install_github("astaaudzi/mizer-rewiring", ref = "temp-model-comp")
#documentation here:
#https://sizespectrum.org/mizer/dev/index.html
library(mizer)
require(tidyverse)
## Loading required package: tidyverse
## Warning: package 'tidyverse' was built under R version 3.6.3
## -- Attaching packages --------------------------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.0
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## Warning: package 'ggplot2' was built under R version 3.6.3
## Warning: package 'tibble' was built under R version 3.6.3
## Warning: package 'tidyr' was built under R version 3.6.3
## Warning: package 'readr' was built under R version 3.6.3
## Warning: package 'purrr' was built under R version 3.6.3
## Warning: package 'dplyr' was built under R version 3.6.3
## Warning: package 'stringr' was built under R version 3.6.3
## Warning: package 'forcats' was built under R version 3.6.3
## -- Conflicts ------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
remotes::install_github("sizespectrum/mizerExperimental")
## Skipping install of 'mizerExperimental' from a github remote, the SHA1 (2e116191) has not changed since last install.
## Use `force = TRUE` to force installation
library(mizerExperimental)
Next, let’s load in and check the data used in the ZooMSS model. TODO: add in some biomass data
#read in group data
groups <-read.csv("data/TestGroups_mizer.csv")
groups$w_min <- 10^groups$w_min
groups$w_inf <- 10^groups$w_inf
groups$w_mat <- 10^groups$w_mat
#groups<-merge(groups,dat,by.x="species",by.y="group",all=T)
# have a look at plot
# plot(groups$w_inf,groups$biomass.tperkm2,xlab="log Maximum Weight [g]", ylab=" log Total Biomass", log="xy",col="blue",pch=16)
# text(groups$w_inf,groups$biomass.tperkm2,labels=groups$species,cex=0.5)
# could plot the paramter allometries here to explore
Next let’s read in the parameters from ZooMSS.
#read groups again for southern ocean model, this time subsetting key groups
groups <-read.csv("data/TestGroups_mizer.csv")
groups$w_min <- 10^groups$w_min #convert from log10 values
groups$w_inf <- 10^groups$w_inf
groups$w_mat <- 10^groups$w_mat
groups$h <- 100000 #maximum feeding level - very large value to reflect the unlimited feeding implicitly allowed in ZooMSS
#groups <- readRDS("data/groups.RDS")[-1,]
#check fails these tests:
groups$w_mat25 >= groups$w_mat #note we don't have w_mat25. Not sure if it's necessary/desirable. Probably not.
## logical(0)
groups$w_mat25 <= groups$w_min
## logical(0)
groups$w_inf <= groups$w_mat25
## logical(0)
#[-1,]
#fix one value
#groups[7,"w_min"]<-4000
# read interaction matrix
# get the interaction matrix - actually I think we can leave this out. Default is all 1s, which is the same as in ZooMSS. Included for completeness; it may be useful in future to keep this in.
theta <- readRDS("data/zoomss_inter.rds")[,-1]
#[-1,-1]
We will pass these parameters to mizer to set up a new multispecies model.
TODO: adjust parameters here.
mf.params <- newMultispeciesParams(species_params=groups,
interaction=NULL, #NULL sets all to 1
no_w = 178, #number of zoo+fish size classes;
kappa = 1e15, #set huge to recreate unbounded feeding on phyto in ZooMSS
min_w_pp = 10^(-14.5), #minimum phyto size
w_pp_cutoff = 10^(-8), #maximum phyto size
r_pp = 20, #Coefficient of the intrinsic resource birth rate, set large for now
n = 0.7, #The allometric growth exponent used in ZooMSS
z0pre = 1, #external mortality
z0exp = 0.3,
resource_dynamics = "resource_constant",
#RDD = constantRDD(species_params = groups) #first go at this
#pred_kernel = ... #probably easiest to just import this/pre-calculate it, once dimensions are worked out
)
## Note: No ks column so calculating from critical feeding level.
## Note: Using z0 = z0pre * w_inf ^ z0exp for missing z0 values.
# # mf.params@species_params$metab[1] <-0.01
# mf.params@species_params$h[] <- 10*mf.params@species_params$ks[]
#
# # mf.params@species_params$gamma[1] <- 1000*(mf.params@species_params$gamma[1])
# mf.params@species_params$alpha[] <- 0.6
# #mf.params@species_params$z0[1] <- 0.01
#
# #higher mean betas
# mf.params@species_params$beta[mf.params@species_params$species == "large.divers"] <-6000
# mf.params@species_params$beta[mf.params@species_params$species == "apex.predators"] <-1000
# # # wider feeding kernels
# mf.params@species_params$sigma[] <- 2
# mf.params@species_params$sigma[mf.params@species_params$species == "baleen.whales"] <- 4
#
#
# # mf.params@species_params$w_mat25[1] <- 0.005*mf.params@species_params$w_inf[1]
# mf.params@species_params$erepro <- ifelse(mf.params@species_params$erepro==1,0.001,mf.params@species_params$erepro)
#
# # mf.params@interaction[] <-0.1
# #mf.params@interaction[,1] <-0
# #mf.params@interaction[1,1] <-1
#
# #mf.params@species_params$interaction_p[2:11] <-0.5
#
# # hand tuning to gte in line with data
# mf.params@species_params$R_max <- 10*mf.params@species_params$R_max
#
#
# # needs to be very different for marine mammals? use different repro assumptions (fixed offspring density per year)?
#
#
#
# # mf.params@species_params$w_mat <- 0.5*mf.params@species_params$w_inf
# # mf.params@species_params$w_mat[mf.params@species_params$species == "flying.birds"] <- 0.75*mf.params@species_params$w_inf[mf.params@species_params$species == "flying.birds"]
# # mf.params@species_params$w_mat[mf.params@species_params$species == "small.divers"] <- 0.75*mf.params@species_params$w_inf[mf.params@species_params$species == "small.divers"]
# # mf.params@species_params$w_mat[mf.params@species_params$species == "medium.divers"] <- 0.75*mf.params@species_params$w_inf[mf.params@species_params$species == "medium.divers"]
# # mf.params@species_params$w_mat[10:12] <- 0.75*mf.params@species_params$w_inf[10:12]
# # setParams(mf.params)
#
# mf.params<- setParams(mf.params,kappa=1e5) # take the new paramewters and change kappa
sim <- project(mf.params, t_max=500,dt = 0.1)
plot(sim)
#plotlyGrowthCurves(sim,species="macrozooplankton")
plotlyFeedingLevel(sim)
# feeding level satioation for some groups, except for the seabirds
# macrozooplankton - they are not growing enough,why?
#tuneParams(mf.params)
plotlyGrowthCurves(sim,percentage = T)
## DLSODA- At current T (=R1), MXSTEP (=I1) steps
## taken on this call before reaching TOUT
## In above message, I1 = 5000
##
## In above message, R1 = 0.127549
##
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
## DLSODA- At current T (=R1), MXSTEP (=I1) steps
## taken on this call before reaching TOUT
## In above message, I1 = 5000
##
## In above message, R1 = 0.00198096
##
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
## DLSODA- At current T (=R1), MXSTEP (=I1) steps
## taken on this call before reaching TOUT
## In above message, I1 = 5000
##
## In above message, R1 = 0.00198096
##
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
plotlySpectra(sim)
—–Old stuff only below this line—–